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Abstract 

This paper introduces a reaction-diffusion-chemotaxis model for bacterial aggre- 
gation patterns on the surface of thin agar plates. It is based on the non-linear 
degenerate cross diffusion model proposed by Kawasaki et al. [Ij and it includes 
a suitable nutrient chemotactic term compatible with such type of diffusion. 
High resolution numerical simulations using Graphic Processing Units (GPUs) 
of the new model are presented, showing that the chemotactic term enhances the 
velocity of propagation of the colony envelope for dense-branching morphologies. 
In addition, the chemotaxis seems to stabilize the formation of branches in the 
soft-agar, low- nutrient regime. An asymptotic estimation predicts the growth 
velocity of the colony envelope as a function of both the nutrient concentration 
and the chemotactic sensitivity. For fixed nutrient concentrations, the growth 
velocity is an increasing function of the chemotactic sensitivity. 

Keywords: Nutrient chemotaxis. Bacillus subtilis, nonlinear degenerate 
diffusion, front propagation. 



1. Introduction 

The spontaneous emergence of complex patterns is ubiquitous in nature. 
When faced against hostile environmental conditions, for example, bacterial 
colonies may exhibit very diverse morphological aggregation patterns. Such 
conditions can be reproduced during in vitro experiments on the surface of thin 
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Figure 1: Depiction of the morphological diagram of the aggregation patterns exhibited by 
colonies of B. subtilis as a function of the concentration nutrient Cn (vertical axis) and the 
solidity of the agar medium expressed as l/Ca (horizontal axis). Both concentrations are 
measured in grams per liter [g/1]. This diagram was redrawn from the one in 1 . 

agar plates by using a very low level of nutrients, a hard surface (high concen- 
tration of agar), or both. Some species (like the bacterium Bacillus subtilis) 
are known to exhibit a great variety of branching patterns depending on the 
nutrient concentration, softness of the medium, and temperature. For example, 
when inoculated on a nutrient-poor solid agar, the patterns show fractal mor- 
phogenesis similar to diffusion-limited aggregation (DLA) (see, e.g., Matsuyama 
and Matsushita [2^, and Ben- Jacob et al. [3 ). When inoculated on semi-solid 
agar (softer medium) the bacterial colonies exhibit a dense-branching morphol- 
ogy (DBM) with a smooth colony envelope propagating outward (cf. Ohgiwari 
et al. [4 ). Finally, when both nutrient and agar's softness further increase, 
the bacteria form simple circular patterns growing homogeneously (cf. Wakita 
et al. [5 ). A detailed account and comparison of the different experimental 
results can be found in Kawasaki et al. [1 (see also Golding et al. [6 and the 
references therein). These experimental observations are usually summarized in 
a morphological diagram like the one shown in Figure [l] The horizontal axis 
measures the medium softness expressed as 1/Ca, where Ca is the agar con- 
centration, and the vertical axis is the concentration of the nutrient Cn- When 
the medium is very hard, the patterns are DLA-like (region A) at low nutrient 
concentration, and round and compact with fractal boundary for high values of 
Cn (region B). DBM-like patterns arise in regions where the nutrient is poor 
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and the agar softness is intermediate (semi-solid), with an envelope smoothly 
rounded and moving outward (region E). In region D, where the agar is softer 
and the nutrient is more abundant the dense branches fuse together forming 
homogeneous colonies with smooth boundaries. Although susceptible of further 
improvements (see, for example, Wu, Zou and Jin [7 ), these diagrams provide a 
good macroscopic description of the colony patterns in terms of simple physical 
incubation conditions. 

In order to explain the experiments, many mathematical models have been 
proposed. One method of modeling is the continuous deterministic reaction- 
diffusion approach in which the bacterial density and the nutrient concentra- 
tion are described with continuous time evolution systems of partial differen- 
tial equations (cf. Murray [8]). Motivated by the experimental observations of 
Ohgiwara et al. 0] for DBM-like transition patterns between regions E and D in 
the morphology diagram (see section 2.1 below), Kawasaki et al. [1 proposed 
in the mid-90's a reaction-diffusion model with non-linear density-dependent 
degenerate (cross) diffusion coefficient, which captures many of the pattern fea- 
tures found experimentally. Thanks to its capability of reproducing the complex 
dense morphology of patterns as well as its rich mathematical structure, this 
model has established itself as a well-known non-linear diffusion system in the 
literature (see, e.g., [8l|9l[T0] and the references therein). 

There are, however, other phenomena which must be taken into consider- 
ation. Bacterial chemotaxis (that is, the process by which bacteria migrate 
toward higher concentrations of certain chemical fields) has attracted signifi- 
cant interest due to its critical role in pattern formation (see [TTJ [121 HSl HI] and 
[15] for a review). Ben- Jacob and co-workers (cf. [El El Hi]) have suggested 
that the chemotactic response of the bacteria is essential to understand some 
of the fundamental features of the observed patterns. In particular, Golding et 
al. [6] discussed the experimental evidence and provided theoretical arguments 
favouring the introduction of chemotaxis towards high concentration of nutri- 
ents into many of the the reaction-diffusion models existing in the literature. 
Although the authors did not propose a new chemotactic term for the particular 
model of Kawasaki and collaborators, they certainly discussed what the relation 
between a generic diffusion term and the appropriate chemotaxis should be (see 
section 



2.2 below). 



Following these ideas, this paper incorporates a suitable nutrient chemotac- 
tic term into the original model by Kawasaki et al.^ and explores the effects 
that such chemical signaling brings upon the aggregation patterns. For that 
purpose, we have conducted high resolution numerical simulations of the new 
system of equations using Graphic Processing Units (CPUs) to perform par- 
allel computations. Moreover, we were able to reproduce the results without 
chemotaxis and to compare them to the new emerging patterns when the chem- 
ical signaling is switched on. One of the main features of the new chemotactic 
term is an increase of the front propagation velocity for the colony envelope. 
Such observation is justified using detailed front asymptotic calculations of the 
speed and the numerical results. We find that, except for a regime of soft agar 
and low nutrient, the chemotaxis does not modify the pattern colony structure 
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significantly, providing evidence of the conjecture of Cohen et al. {16* . 

The remainder of this article is organized as follows. Section [2] gives ac- 
count of the original non-linear reaction-diffusion (non-chemotactic) model and 
introduces the chemotactic term into the equations. The resulting system is 
further normalized in order to reduce the number of parameters. Section [3] con- 
tains the results of our numerical simulations, as well as a comparison of the 
observed patterns to those without chemotaxis. In order to explain the effects 
of the chemotactic terms, section [4] contains an approximation of the speed of 
the envelope front using the ideas of geometric front propagation. Furthermore, 
we provide a numerical estimation of the speed which approximates well the 
theoretical result, and compare it to the speed without chemotaxis. The es- 
timations are supported by some theoretical calculations which can be found 
Appendix A at the end of the paper. Finally, in section [5j we make some 



m 



concluding remarks. 



2. Mathematical modeling 

In this section we introduce the mathematical reaction-diffusion-chemotaxis 
model studied in this paper. We start by providing a detailed explanation of 
the original model in pLj. 

2.1. The model of Kawasaki et al. 

The model proposed by Kawasaki et al. [1] is a reaction-diffusion continu- 
ous system of partial differential equations that takes into account the detailed 
movement of the bacteria. The model has the following general form: 

fit = DnAn - f{n, h) (la) 
ht = V -{D^Vh)^ef{n,h) (lb) 

where n = n(x, t) and h = 6(x, t) represent the concentration of the nutrient 
and the density of the bacterial cells, respectively. We will denote every point 
in the domain as {x^y) G 11 C . Here (1 is a bounded domain and t > 
denotes time. This model assumes that the bacteria and the nutrient diffuse, 
with diffusion coefficients Di^ and D^^ respectively. The coefficient Dn > is 
constant. 

The key feature of this model is the choice of the nonlinear diffusion coef- 
ficient 1^5, which depends on both h and n. Its design was motivated by the 
work of Ohgiwari et al. [4 who investigated the patterns formed by B. subtilis 
on thin agar plates. They discovered that variations in the environmental con- 
ditions as well as in the concentrations of both nutrient and agar medium lead 
to drastic changes in the morphology of the colony patterns. They found two 
distinct type of growing processes. One is performed by immobile cells when 
the growth occurs on hard agar plates. In this case, the outermost part of the 
colonies grows by cell division and cells in the inner part change into a dormant 
state. The other type of growth takes place for intermediate and low values of 
Ca- This means that the growth is produced by the movement of bacterial cells 
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and, due to the softness of the agar, they can move around dynamicahy. The 
outermost front of the colonies is a boundary layer of cells whose movement 
is dull. But cells inside the layer move around actively and sometimes they 
break through the layer and immediately become dull. The bacteria inside the 
colonies are usually inactive. 

In view of such experimental observations, Kawasaki and co-workers conjec- 
tured that the bacteria are immotile when either n or h are low and become 
active as n or 6 increase. Therefore, they proposed Di^ as a non-linear (cross) 
diffusion function of the form 

= crnh, (2) 

where a = <Jo(l + A). The parameter <Jo > is constant and represent the agar 
concentration (its value increases as the concentration of the agar decreases), 
whereas the parameter A captures the random diffusion of the cells as a random 
fluctuation. 

Finally, the term /(n, 6) is the consumption rate of the nutrient by the 
bacteria and ^/(n, h) is the growth rate of cells, being ^ > the conversion rate 
factor. It is assumed that the function / is of Michaelis-Menten type, 

where > and 7 > are considered as constants. In the case when the 
concentration of the nutrient is low, the function ([3| can be approximated by 

/(n, h) = nnh. (4) 

The latter is the form of the kinetic function used by Kawasaki and collaborators 
in their simulations, and the function that we consider in this paper as well. 



2.2. The role of chemotaxis 

Golding et al. [6 argue that the growth of bacterial colonies is an auto or- 
ganization process which involves several mechanisms of chemotactic responses. 
They suggested that there are at least three different types of chemotactic sig- 
nals. One of the signals is the result of nutrient (or food) attractive chemotaxis 
and it is the dominant signal in certain regimes of the nutrient level. Another 
type of signal is a long range repulsive one, which is produced by the starving 
bacteria in the center of the colonies. The last type of signal is secreted by the 
bacteria in the front that are immersed in waste products. These bacteria ask 
for help to metabolize waste products by segregating a chemoattractant. The 
range of this signal is determined by the diffusion coefficient and the rate of 
decomposition. It is, in general, of short range. 

When refering to bacterial colonies, Golding et al. [6 suggested that nutrient 
chemotaxis is the mechanism responsible of the increment in the growth speed 
and the control (or even the decrease) of the fractal dimension of the arising 
patterns. Such process should permit to increase the propagation of the front. 
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In view of these observations, we introduce a nutrient chemotactic term into 
system ([T]). For that purpose, we need to define the chemotactic fiux Jc which, 
in general, is written as 

Jc = mx{n)Vn. (5) 

Here x = x(^) is the chemotactic sensitivity function and ^ = ^(6) is the bac- 
terial response to the nutrient gradient. The chemotactic signals are detected 
by dedicated membrane receptors which perform a binding between ligands (at- 
tractors) and specific binding proteins [18]. Bacteria sense temporal gradients, 
that is, they compare sequential registers of their occupied chemoreceptors [19 . 
This complexity is built into the model equations through a signal-dependent 
chemotactic sensitivity function. One of the most commonly used forms is the 
"receptor law" proposed by Lapidus and Schiller [20 , 

where xo > is a constant measuring the strength of the chemotaxis and 
Kd > is the receptor-ligand binding dissociation constant, which has nutrient 
concentration units, and represents the nutrient level needed for half of receptor 
to be occupied. This receptor sensitivity law has been derived and applied in 
several models for chemotaxis (cf. [21] [22] [23j ) . Its design was motivated by the 
experimental observation of Mesibov, Ordal and Adler that for very high 
concentrations of nutrient the chemotactic response vanishes due to saturation 
of the receptors. This feature was not previously captured by the standard 
Keller-Segel chemotactic fiux function (see 23)- The functional form of (|6| can 
be understood in terms of a simple model for receptor-signal binding through the 
law mass action (see [25l|24]). It is to be noted that the dissociation constant Kd 
for the attractor-receptor molecular interaction has a unique value that depends 
on the bacterial strain, nutrient type and other in vitro conditions, and has to 
be determined experimentally (see Goldman and Ordal [26 for values of Kd 
on various substrates of B. Suhtilis). Moreover, it can be interpreted as the 
nutrient level where the chemotactic sensitivity x(p) is maximum. Thus, we 
shall use Kd as the normalization factor for the nutrient density n. 

The bacterial response function ^ = ^(6) is positive for attractive chemotaxis 
and negative for a repulsive one. In the authors suggest that, if the bacteria 
move within a liquid and the colony density is low, then the bacterial response 
should be proportional to the bacteria concentration times the diffusion, that 
is, 

m\^hDh- (7) 

Notice that, on account of ([7|), if the diffusion coefficient is constant then the 
bacterial response function is proportional to 6, recovering in this case the clas- 
sical Keller-Segel chemotactic ffux ^T1[T2]. Therefore, the authors suggest that, 
for general non-linear diffusion functions the expression of the chemotactic 
fiux Jc must be modified according to ([7| in order to incorporate the density 
dependence of the bacterial movement. For example, for the non-linear diffusion 
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considered by Kitsunezaki [27], namely Di^ = Dob^ where I^o is constant and 
m > is an integer, they proposed ^(6) = bDi^ = DQb^~^^ as the appropriate 
response function. 

2.3. Incorporation of the nutrient chemotactic term 

We now introduce a suitable chemotactic term into model ([T]) expressing the 
chemical signaling of the nutrient concentration n. In view of the form of the 
non-linear diffusion ([2| and taking into account ([7|, we propose the following 
(attractive) chemotactic flux 

Jc = crnb^x{n)Vn, (8) 

where x(^) is the receptor law ( [g]). Substracting the divergence of such flux to 
the bacterial density equation in (ll| we obtain the following reaction-diffusion- 
chemotaxis model: 

rit = DnAn — nnb (9a) 

bt = V • {(jnbVb) + Onnb -V ■ ( an6 ^ J^^^ , Vn^ , (9b) 

V + ny J 

for (x, ?/) G 1^, t > 0. The system is subject to initial conditions modeling an 
uniformly distributed initial constant concentration of nutrient and the initial 
inoculum of the bacteria. They have the form 

n(x, y, 0) = no, h{x, y, 0) = 6o(x, y). (10) 



Here no is a positive constant and b^ is a given function. The system ([9|) is 
further endowed with no- flux boundary conditions, 

Vb-v = {), Vn-i/ = 0, {x,y)edVt, (11) 

where G M^, = 1, is the outer unit normal vector at the boundary of the 
domain. 

System of equations ([9| adds the nutrient chemotactic signaling to model ([T]). 
The flux (|8| can be interpreted as the cross chemotactic function corresponding 
to the nonlinear diffusion ([2|. 

2.4- Normalization 

In order to reduce the number of parameters in system (|9| we introduce the 
following rescaling of the variables: 
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Upon substitution into ^ and dropping the tilded notation for simplicity we 
obtain the fohowing dimensionless system of equations 



nt 
bt 



An — nb 
V • {(jnhVh) 



nb - xoV • 



an 



62 



{i + ny 



with a = (Jo(l 



A), and with initial conditions 

n{x,y,0) = ho/Kd. = no, 

b{x, y, 0) = bo{x, y, 0)/{eKd) = bo{x, y, 0), 



(12a) 
(12b) 

(13) 



for (x^y) G ft. To complete the system we impose no-flux boundary conditions 
on the rescaled variables 



(14) 



The only remaining free parameters are ctq > 0, measuring the hardness of 
the agar medium (large ctq for low agar concentration); no > 0, measuring the 
relative initial nutrient concentration; and xo ^ 0? which measures the intensity 
of the chemotaxis. 

In the sequel we compute numerically the solutions to system (12) under 
conditions ( 13 ) and ( 14 ). We are particularly interested in the interplay between 
the different parameter values expressing agar concentration, initial nutrient and 
the chemotactic signal, and how they affect the colony patterns. 



3. Numerical simulations 



We performed several numerical simulations of the system (12) using an 
explicit second order Runge-Kutta and finite differences numerical scheme and 
the NVIDIA CUD A libraries for the C++ programming language. The libraries 
are written to run on CPUs, and the CUDA programming language is designed 
under the simple program/multiple data (SPMD) programming model (cf. [28 ). 
More specifically, the GPU architecture is very well-suited to address problems 
that can be expressed as data-parallel computations, that is, the same program 
(finite difference calculations) is executed on many data elements (grid points) 
simultaneously (see [29j). Hence, we overcome the inconvenience of small time 
steps required by explicit numerical schemes for solving stiff equations, and were 
able to accomplish thousands of iterations in a couple of hours. Our numerical 
experiments were performed on a NVIDIA Tesla® C2070 graphics card with 
448 CUDA cores. 

The calculations were performed on a two dimensional square domain of 
side L = 680. The selected mesh for all our computations is a square grid of 
N = 2048 points per side (see discussion on section [3^ below) . Therefore, the 
grid width that we used was Ax = Ay = L/N ^ 0.3320. We considered two 
values for the time step, namely. At = 0.0231 and At = 0.0110, for which the 
method showed to be stable. 
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Description 


Symbol 


Values 


Initial condition for the nutrient 


no 


1.07, 0.71 


Agar softness 




1.0, 4.0 


Chemotactic signal strength 


Xo 


0, 2.5, 5.0, 7.5, 10.0 



Table 1: Free parameter values used during the numerical simulations. 



Following [1 and for the sake of comparison, the initial distribution for the 
bacteria was taken as 

bo{x,y,0) = bMe-^-'+y"y'^-^'^ (15) 

where 6m = 0.71 is the maximum density at the center of the domain. Recall 
that the initial condition for the nutrient is a fixed constant, for which we take 
the values no = 0.71 and 1.07. 

Similarly as the computations in [1 , and in order to incorporate the fact 
that bacteria move according to a biased random walk, on every grid point and 
at each time step the diffusion coefficient for the bacteria a is perturbed from 
the mean ao through a random variable A taken from a triangular distribution 
with support on [—1, 1], that is, <j = (Jo(l + A). Remember that (Jq represents 
the softness of the agar medium. 

Regarding the chemotactic term, the coefficient xo represents the strength 
of the signal. The simulations were performed using several values of this pa- 
rameter including the value xo = 0, that is, in the absence of chemotaxis. 
Thus, we also computed solutions to the original system ([T]) as well (in its non- 
dimensional form) also for the sake of comparison. For a complete summary of 
the free parameter values used during the numerical experiments, see table [T] 



3.1. Results 



We now present the numerical simulations of system (12) subject to initial 
conditions (15) and n(x,7/, 0) = no, as well as to boundary conditions (14). We 
divide the exposition of our numerical results in terms of the agar's hardness. 
The latter is expressed by the parameter value do. In our simulations, we 
considered the values ao = 1 and do = 4, and the observed patterns correspond 
to regions E and D in the morphology diagram (see the enclosed region inside 
a rectangle in figure [T]). For do = 1 we took the time step At = 0.0231 and for 
(Jo = 4 we took it as At = 0.0110. The computed solutions are represented in 
figures [5] through [9| Each figure contains two columns associated to two different 
values of the parameter xo^ and each row is associated to a different time step. 
Such configuration allows the reader to compare the colony patterns for different 
chemotactic intensities. In what follows we highlight that the main effect of the 
chemotaxis is the increment on the growth speed of the outer envelope of the 
bacterial colony. 
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3.1.1. Semi-solid agar 

We begin with the case of semi-sohd agar, for which ctq = 1. This means 
that the concentration of agar is intermediate and the diffusion of the bacteria is 
hmited. The observed patterns correspond to the DBM-hke branching structure 
in region E of the morphology diagram. The considered values for the initial 
nutrient level were no = 0.71 and 1.07. 

Figure [5] shows the numerical calculations for no = 0.71. On the left column 
we find the solutions to the system in the absence of chemotaxis, that is, when 
Xo = 0. The pictures for this case are comparable to those obtained by Kawasaki 
et al. (see figure 3(b) in [1 ). On the right column we find the case where the 
chemotactic sensitivity is set as xo = 5. The morphology shown in both cases 
is very similar (if we look carefully we can see that the number of branches is 
consistent). There is, however, a significant increase in the growth velocity of 
the outer envelope of the branches when the chemotactic signal is switched on, 
as the reader may observe. 

The simulations for no = 1.07 are depicted in figure [6j Once again the 
left column corresponds to the value Xo = 0. The patterns are defined by fine 
radially oriented branches with round envelope. At the center of the colony 
there is the highest concentration of bacteria. Actually, the peak concentration 
is the remainder of the initial inoculation, because the bacteria at early times 
deplete the nutrient in site of the initial condition, but the cells in the inner 
part of the colony become immotile. Meanwhile, the cells at the periphery of 
the colony move actively searching for higher concentrations of nutrient; this is 
what drives the formation of branches. Hence, bacteria at the front of the colony 
account for the instability at the interface. This is the main feature of the non- 
linear cross diffusion model. In the right column of figure |6] we find the patterns 
in the presence of a chemotactic signal (xo = 5.0). We observe that at early 
stages (see figure |6(b)| ) the chemotaxis provides a significant outward drift to 
the cellular movement as discussed in [6 , causing a more ramified pattern with 
round envelope and a substantial enhancement of growth speed. Comparing 



figures 6(g) and 6(h) we notice that the chemotactic signal produced a pattern 



two times bigger (the envelope is moving faster). 
3.1.2. Soft agar 

In our next numerical experiment, we set the value ao = 4, causing the 
diffusion of the bacteria to be increased by lowering the concentration of the 
agar. At sufficient high levels of nutrient the colony patterns are restricted to 
region D in the morphology diagram, whereas at lower levels of no patterns 
in region E emerge. In order to explore the morphology transitions between 
regions D and E, we first set no = 0.71 as the initial nutrient and consider 
various values for xo, because the chemotactic signal gives rise to a outward 
drift to the cellular movement. The computations are depicted in figures [7| 
and [S] In this experiment, the comparison is made through the two figures, 
for which the initial nutrient is fixed and the chemotactic sensitivity takes the 
values Xo = 0,2.5,5.0 and 7.5. In figure [7) the left column shows the evolution 
of the system without chemotaxis. At early times the bacteria evolve as a solid 
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colony (figure |7(a)[ ) , but for the consequent iterations they develop a pattern 
of radially oriented fjords. On the right column of figure [7| we show the results 
with chemotactic strength xo = 2.5. As expected, the chemotactic term makes 
the colony grow faster; compare, for example, figures |7(f)| and 7(g) Observe 
that the pictures in the two columns of figure [t] (with chemotactic sensitivities 
Xo = and xo = 2.5) exhibit similar morphologies, being the envelope velocity 
the only apparent difference between them. 

For higher values of xo, however, the growth speed is not only greatly in- 
creased, but there is also a change in the morphologies as well. We may appre- 
ciate, for instance, that for xo = 7.5 the colony has almost doubled its size at 
time step t = 324.0 (see figures 7(g) and |8(h)| ). In this case the increment on the 
growth speed is accompanied by a notable morphological change. Presumably, 
at higher values of xo the propagation speed of the front is higher than the 
diffusion of the nutrient preventing the fjords from forming. Therefore, the thin 



branches start fusing (figure 8(g) ) until the pattern shows no more branches 
and becomes an homogeneous disk (figure |8(h)[ ) . These observations suggest 
that, in the transition between regions E and D in the morphology diagram, 
the supressing effect of the chemotactic term on the onset of instability (which 
causes the formation of branches; cf. [30 ) is more evident than in other regions 
of the diagram. An increase of the chemotactic intensity prevents the formation 
of further branches and results in an homogeneization of the colony morphology, 
accompanied with the expected increase in the velocity of the overall colony. 

Finally, in figure |9] we present the results of our numerical simulations for the 
value of nutrient no = 1.07 and chemotaxis sensitivities Xo = and xo = 7.5. In 
this case no significant morphological changes are observed. In the left column 
of figure [9] the pattern shown is an homogeneous disk with a peak of bacterial 
concentration at the center. In the right column (when xo = 7.5), as expected, 
it is observed that the chemotactic signal increases the growth speed and the 
pattern is also a slightly rugged solid disk, but with a small depression at the 
center. The latter happens because the softness of the agar allows the colony 
to spread out before the nutrient is depleted. In this case (which correspons to 
region D in the diagram), we notice the increment in the outer colony velocity 
with no significant change in morphology. 

3.2. A note on the grid size 

It has been suggested that numerical simulations might depend on the grid 
size [31 . To circumvent this possibility, we performed simulations with the 
following grid sizes as weh: N = 1024,2048,4096 and 8192. The results for the 
parameter values ctq = 1, no = 0.71 and Xo = 5 can be observed in figure [2] for 
each grid size. Notice that in the case of = 1024 points per side, the solutions 
exhibit a cross shaped pattern with a rhomboid envelop (figure 2(a)). When 
the grid size is increased to A" = 2048 there is indeed a significant change in 
morphology, showing a denser branching pattern with a round envelope (figure 
|2(b)[ ). This pattern does not change significantly, however, when the size grid 
is increased to A" = 409 6 and A = 81 92 (fig ures |2(c)| and [2(d)| . The reader 
may notice that figures 2(b), 2(c) and 2(d) show very similar patterns. We 
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conclude that a resolution of N = 2048 lies above the threshold needed for a 
reliable display of the solutions of the system. Even though our computations 
were performed with parallelization techniques on a GPU, the simulations with 
N = 8192 involved a large computational time (there are more than 64 million 
grid points). Therefore, we adopted N = 2048 as the resolution for all the 



simulations of system (12) in this paper. 




(a) N = 1024 



(b) N = 2048 




(c) = 4096 



(d) = 8192 



Figure 2: Simulations of system ( |12[ ), with parameter values ctq = l,no = 0.71 and xo = 5, 
for different values of the grid size A^. 



4. Propagation velocity of the front 

In this section we discuss the effects of the chemotactic term on the envelope 
propagation speed by means of an asymptotic geometric front approximation. 
Moreover, we provide numerical estimations of the speed by simulating a one 
dimensional version of system (12). We compare these velocity calculations to 
those without che mot axis. 
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4.1. Approximation of the speed 

For simplicity, let us consider the chemotactic sensitivity function as a con- 
stant, i.e., X = Xo- Therefore, the system of equations takes the form 

fit = An — nb (16a) 
ht = V ' {(jnhVh) ^ nh - xo^rV • {n}?Vn) , (16b) 

In order to approximate the velocity of the smooth envelope enclosing the 
branching patterns, we explore the following reduction proposed by Kawasaki et 
al. [Tj . They made a further approximation in order to obtain a scalar equation 
for the bacterial concentration h. They observed that in the absence of diffu- 
sion the total mass is conserved. In our setting, if we ignore both diffusion and 
chemotaxis then we can add the equations and integrate them in time to obtain 

n^h = C (17) 

where C is a constant, which will be taken as C = no- Thus, n = no — b and we 
obtain the logistic equation 

bt = nob f 1 - 

V ^0 

It must be observed that in [l] the term aonb is substituted by aono^, and 
the term nb is substituted by no6(l — 6/no), in the approximation of the system 



by a scalar equation for b. Instead, we are substituting (17) into (16b) to obtain 
a scalar equation for b. This not only represents a more accurate approximation 
but it is also consistent with the numerical computations which, at first order. 



show that the conservation- like equation (17) is valid away from the front. 



Substitution of ( |T7| ) into ( |T6| ) yields the scalar equation 

6t = V • (ao(no - b)b\/b) + (no - b)b - aoXoV • ((no - 6)6'V(no - 6)), (18) 

which is an approximated scalar equation for b with a Fisher-type reaction 
term. Here we have approximated <Jo(l + A) ^ do, because in this part of the 
study we are not interested in considering the random motion of bacteria. A 
rearrangement of the terms in equation ([T8| allows us to write the two divergence 
terms together. This yields 

bt=aoV-{D{b,xo)Vb)+g{b), (19) 

where 

D{b,xo)=nob(^l- ^^{l + Xob), (20) 

may be interpreted as an effective non-linear diffusion coefficient, which involves 
both the contributions of diffusion and chemotaxis. The reaction term is of 
Fisher-KPP type, 

g{b) = nob (l - • (21) 
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We describe the motion of the front in local curvilinear coordinates with 
components ({x^y^t) and r{x,y,t), the normal and tangent unit vectors to the 
front, respectively. According to custom, ( is the normal component pointing 
inside the colony front, so that —Q is the outer normal speed. We assume that 
the dependence of b with respect to the tangent component is negligible and we 
approximate 

b{x,y,t)^b{C{x,y,t)). (22) 



After substituting b into equation (19), and dropping the bar for notation con- 
venience, we get the following ordinary differential equation for b 

- cb' = ao{D{b)by + b{no - b). (23) 

where 

c=-Ct+aoDib)AC (24) 

and —(t is the velocity of the envelope front that we want to estimate. 

The approximated front equati on ([23| ) in the normal direction resembles a 
one-dimensional equation of form ( A.6| ) (see Appendix A). By the results of 



Malaguti and Marcelli [32j the one-dimensional velocity for each value of xo > 
can be estimated in terms of Xo, ^o, and the functions D and g. We note that we 
are interested in comparing the velocity between the cases xo = and Xo > 0- 
The aforementioned result shows that one-dimensional traveling wave solutions 



to equations of the form (A.6) must travel with velocity c(xo) bounded by 



c(xo) ^ ^^(Xo), where the threshold velocity satisfies the bound: 



< c*(xo) < c(no,Xo) := 2W sup (25) 

V beiO,no] ^ 



with D and g given by (20) and (21 ). Upon substitution, we are able to compute 
the bound 



b^' 



c(no,Xo) = 2noV^\/ max 6 1 (l + xo^)- (26) 

y 6G[0,no] V ^0/ 

The function ?/;(6) = 6(1 — 6/no)^(l + Xo^) is non-negative for < 6 < no 
and each Xo ^ 0; moreover, it has a global maximum ilJmax = where 




2 V 2 4X0 

and with < 6* < no for all values of xo > and any fixed value of no > 0. 
Thus, 

c(no, Xo) = 2noV^\/V<M- (28) 
In the absence of chemotaxis (xo = 0) the function il){b) reduces to il){b) = 
b{l — b/no)'^ and has a global maximum ?/^(no/3) = 4no/27. This means that 



c(no, 0) = 2noV^o\l ^ = ^^o^'^'^'- (29) 
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Equations (28) and (29) are the speed thresholds for the cases with and without 



chemotaxis, respectively. These equations predict the speed of a sharp front, and 
will be used when comparing the numerical estimates of the velocity (see section 
|4.2|below). In our setting, we are considering that the growth of the colony is in 



the radial direction (see equation (22)), so we only consider the normal velocity 
s = —(t- Furthermore, we take into account the cases where the colony exhibits 
an outer growth, this means that the local curvature is positive there 



-AC > 0. 



(30) 



Returning to equations (23) and (24) and comparing the one dimensional 
velocity of a front for equation (A.6) to expression (24), we obtain 

s = -Ct> c*(xo) + croD{b)K > 0. (31) 

It follows from equation ([30| and Proposition [l] in [Appendix A[ that 

s = -Ct> c*(xo) + (ToD{b)f^ > c*(0) + aoD{b)f^ > c*(0). (32) 

Last inequality implies that the normal velocity s is greater than the ve- 
locity for the sharp front when there is no chemotaxis (xo = 0)^ and that it 
is increased/decreased by a term proportional to the local curvature and the 
chemotactic strength xo > 0- In the outer front where diffusion is degenerate, 
this last curvature term is negligible (as i) = in a vicinity of the envelope 
front) so that we can approximate the speed by the one-dimensional sharp front 
velocity, c*(xo) ^ regardless of the curvature sign. 

4.2. Numerical computation of the velocity 



As we have stated in (22), the tangential direction is negligible in describing 



the front movement. Then, we shall estimate the velocity of the envelope front 



solving numerically a one-dimensional version of system (16), that is of the 
following form. 



dn 

'm 

db 
dt 



dn 
dx 

a^nb 



nb 



db 

dx 



Xocront 



; dn 
dx 



nb. 



(33a) 
(33b) 



The simulations of system (33) were performed on the spacial domain [0,20], 
and the initial conditions were take as 



n(x,0) = no, 6(x,0) = bMe 



(34) 



where 6m = 0.71. Again, as in the two-dimensional simulations, no measures the 
initial level of nutrient, xo measures the chemotaxis intensity and ao measures 
the agar hardness. We made several simulations of system (33) for various 
values of no = {0.35, 0.71, 1.07, 2, 3, 4}, xo = {O7 17 2.5, 5, 10}. In all simulations 
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we fixed ctq = 1. Figure |3] depicts ttie numerical solutions of the bacterial density 
of b at different time steps and different values of no and xo- It should be noted 
that solutions of b initially do not exhibit a traveling wave profile, nevertheless, 
eventually the solution resembles a traveling wave. Therefore, the traveling wave 
estimations were computed at later times by selecting a point in the vicinity of 
0.5. These are the points shown in figure [Sj In figure |4] we present a comparison 
of the propagation velocity of the front between the numerical estimations and 



the velocity threshold given by equations (|28[) for xo > or (29) for xo = 0, as 



functions of the initial nutrient no- For the case without chemotaxis, we found 



sharper results (see figure 4(a) ) than those in the original work by Kawasaki 



et a/., because we are considering a more accurate velocity threshold. 



explained in section 4.1 As expected from equation (32), when a chemotactic 
signal is present the velocity of the front is greatly increased. We can see that 
the numerical velocity fits the velocity threshold. Therefore, figure [4] highlights 
two facts: first, numerical estimations of the propagation speed of the one- 
dimensional system are in good accordance with the theoretical speed thresholds 
that are defined for the scalar equation for b (see equation (Il9|); and second. 



system (16). 



the conservation-like equation (17) is a good approximation of the solutions of 



5. Conclusions 

In this work we have studied the effects that nutient chemotaxis induces 
into bacterial population patterns modeled by a reaction-diffusion system orig- 
inally proposed by Kawasaki et al. [1 . Motivated by a discussion provided by 
Golding et al. [6|, we introduce into system ([T]) a chemotactic term compatible 
with the non-linear cross diffusion. After a suitable rescaling of variables, we 
reduce the number of free parameters of system ([9|, and we only considered ctq, 
no and xo- We solved numerically system ([l2|, for various values of the free 
parameters, including the case without chemotaxis (xo = 0). This means that 
we reproduced the numerical results obtained by Kawasaki et al. in [1 . Our 



numerical simulations of system ( 12 ) show that the chemotactic term provides 
an outward drift to the bacterial movement and therefore enhances the growth 
velocity of the bacterial patterns. 

In order to approximate the velocity of the envelope front, under suitable 
assumptions, we derived a scalar equation of the bacterial density, see equation 



Jl9|). Then, applying a result from Malaguti and Marcelli [32 , see Appendix 
|A[ we proved that the normal velocity s is greater than the velocity of the sharp 
front when there is no chemotaxis (xo = 0)- Afterwards, by simulating a one- 
dimensional version of system ([l6|, we calculated numerical estimations of the 
propagation velocity of the front. We compared the theoretical speed predic- 



tions, defined in equations (28) and ([29]), versus the estimations and we found 
that the propagation velocity of the envelope front is greatly increased when a 
chemotactic signal is present (see figure |4|. For the case without chemotaxis we 



found results comparable with those found by Kawasaki et al. (see figure 4(a) ) 



although we must stress the fact that we are using a more accurate velocity 
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(a) no = 0.71, xo = 10 



(b) no =2,xo = 2.5 




(c) no = 3, xo ^ 



(d) no = 4, xo = 1 



Figure 3: Numerical solution of the bacterial density I 
and different values of no and xo- 



from the system (|33|) at different times 



threshold than the one used in [1 . Furthermore, as figure [4] depicts, the nu- 
merical estimations of the propagation speed of the one-dimensional system are 
in good accordance with the theoretical speed thresholds that are defined for 



the scalar equation for b (see equation (19)), implying that the conservation-like 
equation (17) is a good approximation of the solutions of system (16). 



Finally, we discuss the relation of our results to the analysis of Arouh and 
Levine [30], who studied the Kessler-Levine reaction-diffusion equations (cf. 
[33] ) and explored the addition of a nutrient chemotactic term. They found 
that the chemotaxis suppresses the onset of instability causing the formation 
of branching patterns in the colony envelope as it propagates outward. This 
suppression of instability is refiected in a decrease of the "growth rate" of the 
perturbation of the envelope when the chemotaxis is switched on. The au- 
thors conjecture that this behavior can be extrapolated to other (more realistic) 
chemotactic models like the one introduced in this paper. It is to be pointed 
out, however, that our findings do not contradict the analysis of Arouh and 
Levine inasmuch as our numerical simulations and analytic computations ex- 
hibit an increase in the normal velocity of the envelope outward due to the 
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Figure 4: Comparison of the propagation speeds of the front for the one-dimensional system 
( |33| ) as a function of the initial nutrient concentration for different values of the chemotactic 
sensitivity xo- The asterisks depict the numerical estimations of the velocity. Th e solid line 
is the velocity threshold given by equation \29\ in the absence of chemotaxis ( case|4(a)| >, and 
by equation ( |28| when chemotaxis is switched on (cases |4(b)[ |4(c)[ |4(d)| and |4(e)| , both as 
functions of the initial nutrient concentration no- Here a = 1. 



chemotaxis, whereas the growth rate that Arouh and Levine refer to is the tem- 
poral eigenvalue associated to an eigenfunction (perturbation) of the linearized 
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operator around the envelope front. Actually, a simple calculation using energy 
estimates (not included here) shows that in the case of our model, the spectral 
bounds for the eigenvalues of the linearized operator around the envelope front 
indeed decrease as functions of the chemotactic sensitivity. Our results and the 
observations of Arouh and Levine combined suggest that, although the speed of 
the envelope front increases with the chemotaxis (as it is shown in this paper), 
the growth rate of linear perturbations of the envelope is pushed to the stable 
(negative) region of the spectrum, delaying the formation of fingering unstable 
patterns (as it is shown in [30]). 



Appendix A. Bounds for the velocity 

In this section we summarize the theoretical results of Malaguti and Marcelli 
[32] that allow us to justify why the speed of the envelope front is greater when 
there is a chemotactic signal. Theorem 9 in [32^ states that for one- dimensional 
reaction-diffusion equations of the form 

ut = {D{u)u:,):, + g{u) (A.l) 

with D G C^([0, 1]), ^ G C^([0, 1]), where the diffusion degenerates at 1/ = and 
= 1, namely, 

D{0) = D{1) = 0, D{u) > for ah < < 1, (A.2) 

and the rate of growth of the population is of Fisher-KPP type, 

g{0) = g{l) = 0, g{u) > for all < u < 1, (A.3) 

then there exist a traveling wave solution (unique up to translation) to equation 
(A.l), provided that the velocity c satisfies < < c, where is a threshold 
velocity satisfying the bound 



0<c.<2jsup (A.4) 

ysG(o,i) s 

The case c = corresponds to a sharp front. They prove the theorem by 
showing that the existence of such a traveling wave is equivalent to the existence 
of a solution z = z{u) of the boundary- value problem 

dz D{u)g{u) 

du z ' (A-5) 

z{Q+)=z{l-) = Q, 

for the the same c, satisfying z{u) < of all G (0, 1). 

The result can be easily extrapolated to a one-dimensional equation of the 
form 

bt = {b{b,xM.+9{b), (A.6) 



19 



on [0,no] and where D{b,xo) and g{b) are defined in equations (20) and (21), 
respectively. The existence of a travehng front solution to the equation with 
Xo > is related to the solution to the boundary value problem 

dz^ D{b, xo)ff(6) 

db z ' (A.7) 

^(0+) = z(no-) = 0. 

Likewise the traveling wave without chemotaxis is hnked to the solutions to the 
problem 

dz _ D{b,0)g{b) 

db''"" z ' (A.8) 

z(0+) = z(no) = 0. 

In order to compare the speed of the envelope front in both cases, we consider 
the following sets 



= {c > : problem (A.7) has a non-positive solution}, (A. 9) 



^0 = {c > : problem (A.8) has a non-positive solution}, (A. 10) 

that define the range of values of c for which there exists a traveling wave 
solution for equation ([l9| (see Theorem 9 in [32 ). It should be noted that the 
infimum of the set, in either case, corresponds to the threshold speed c*. The 
following proposition compare the two velocity thresholds (with and without 
chemotaxis). 

Proposition 1. For any fixed value of uq > and for all values of ^0 there 
hold, 

('^) c(xo) > c(0), and 
(ii) c*(xo) > c*(0). 

Proof The proof of (i) is straightforward: for all values of xo ^ and for each 
b G (0,no) there holds 

^{b, xo) = b{l - b/no)\l + xo^) > ^{b, 0) = b{l - b/n^)\ 

Therefore max5^[o,no] V^(^,Xo) > max5^[o,no] V^(^,0), yielding (i). 

To show (ii), let us suppose that c > is a velocity such that c G A^^ for 
some Xo > 0- This means that the problem (A.7) has a non-positive solution 
C{b). Therefore, since D{b^xo) > D{b^O) for all b G (0,no) inasmuch as xo > 0, 
we have that 

dC^_ D{b,xo) D{b,0) 

db ' m ' m ■ 



By Lemma 8 in [32^, there exists a negative solution z = z{b) to problem (A.8), 
with the same constant c. This shows that c G Aq. Therefore, A^^ C Aq for 
all Xo ^ and consequently c*(xo) = i^f ^xo — <^*(0) = inf Aq. This yields 

(ii). □ 
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(a) xo = 0,t = 694.5 



(b) xo = 5, f 



= 694.5 




(a) xo = 0,t = 115.7 (b) xo = 5, t = 115.7 




(c) xo = 0,t = 231.5 (d) xo=^,t = 231.5 




(e) xo = 0,t = 370.4 (f) xo = ^, t = 370.4 




(g) xo = 0,t = 509.3 (h) xo = 5, t = 509.3 



Figure 6: Colony growth as a result of simulations of system ( |12| , taking ctq = 1.0, no = 1.07, 
with chemotaxis sensitivity Xo = (no chemotaxis - left), and xo = 5.0 (right), for different 
values of t. 
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(a) xo = 0, t = 92.6 



(b) xo = 2.5, t = 92.6 
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